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We study a system of TV ^ 1 degrees of freedom coupled via a 
smooth homogeneous Gaussian vector field with both gradient and 
divergence-free components. In the absence of coupling, the sys¬ 
tem is exponentially relaxing to an equilibrium with rate /t. We 
show that, while increasing the ratio of the coupling strength to the 
relaxation rate, the system experiences an abrupt transition from 
a topologically trivial phase portrait with a single equilibrium into 
a topologically non-trivial regime characterised by an exponential 
number of equilibria, the vast majority of which are expected to be 
unstable. It is suggested that this picture provides a global view on 
the nature of the May-Wigner instability transition originally discov¬ 
ered by local linear stability analysis. 


complex systems | equilibrium | model ecosystems | random matrices | 


W ill diversity make a food chain more or less stable? 

The prevailing view in the mid-twentieth century was 
that diverse ecosystems have greater resilience to recover from 
events displacing the system from equilibrium and hence are 
more stable. This ‘ecological intuition’ was challenged by 
Robert May in 1972 [1]. At that time, computer simulations 
suggested that large complex systems assembled at random 
might become unstable as the system complexity increases 
[2]. May’s 1972 paper complemented that work with an ana¬ 
lytic investigation of the neighbourhood stability of a model 
ecosystem whereby N species at equilibrium are subject to 
random interactions. 

The time evolution of large complex systems, of which 
model ecosystems is one example, is often described within the 
general mathematical framework of coupled first-order nonlin¬ 
ear ordinary differential equations (ODEs). In the context of 
generic systems, the Hartman-Grobner theorem then asserts 
that the neighbourhood stability of a typical equilibrium can 
be studied by replacing the non-linear interaction functions 
near the equilibrium with their linear approximations. It is 
along these lines that May suggested to look at the linear 
model 


dyj 

dt 


N 

-m+ '^Jjk Vk, j = i,...,N, 

k=l 


[ 1 ] 


to study the stability of large complex systems. Here J = 
(Jjk) is the coupling matrix and fi > 0. In the absence of in¬ 
teractions, i.e., when all Jjk = 0, system Q is self-regulating: 
if disturbed from the equilibrium yi = y 2 = ■ ■ ■ = yN = 0 it 
returns back with some characteristic relaxation time set by 
fj,. In an ecological context yj (t) is interpreted as the variation 
about the equilibrium value, yj = 0, in the population density 
of species j at time t. The element Jjk of the coupling ma¬ 
trix J, which is known as the community matrix in ecology, 
measures the per capita effect of species k on species j at the 
presumed equilibrium. Generically, the community matrix is 
asymmetric, Jjk 7 ^ Jkj- 

For complex multi-species systems information about the 
interaction between species is rarely available at the level of 
detail sufficient for the exact computation of the community 
matrix and a subsequent stability analysis. Instead, May con¬ 
sidered an ensemble of community matrices J assembled at 
random, whereby the matrix elements Jjk are sampled from 


a probability distribution with zero mean and a prescribed 
variance . This is similar to the approach taken by Wigner 
in his description of statistics of energy levels of heavy nuclei 
via eigenvalues of random matrices, which proved to be very 
fruitful [ 3 ]. A detailed review of May’s model in the light of 
recent advances in random matrix theory can be found in [3]. 

The linear system Q is stable if, and only if, all the eigen¬ 
values of J have real parts less than y. Invoking Wigner’s 
arguments for studying eigenvalues of large random matrices, 
May claimed that for large N the largest real part of the eigen¬ 
values of J is typically ay/N. Obviously, the model’s stability 
is then controlled by the ratio m = y,/{ayfN). For N large, 
system Q will almost certainly be stable if m > 1 and unsta¬ 
ble if m < 1, with a sharp transition between the two types of 
behaviour with changing either y,, a or N. In particular, for 
fixed /r, a. system 0 will almost certainly become unstable 
for N sufficiently large. 

Despite the simplistic character of May’s model , his pi¬ 
oneering work gave rise to a long standing ‘stability versus di¬ 
versity’ debate, which is not fully settled yet n El [71 mis], and 
played a fundamental role in theoretical ecology by prompting 
ecologists to think about special features of real multi-species 
ecosystems that help such systems to remain stable. Vari¬ 
ations of May’s model are still being discussed nowadays in 
the context of neighbourhood stability, see [1] and references 
therein. 

One obvious limitation of the neighbourhood stability 
analysis is that it gives no insight into the model behaviour 
beyond the instability threshold. Hence May’s model has only 
limited bearing on the dynamics of populations operating out- 
of-equilibrium. An instability does not necessarily imply lack 
of persistence: populations could coexist thanks to limit cy¬ 
cles or chaotic attractors, which typically originate from un¬ 
stable equilibrium points. Important questions posing serious 
challenge then relate to classification of equilibria by stabil¬ 
ity, studying their basins of attraction, and other features of 
global dynamics |4]. Over the last years, as the computing 
power grew, non-linear models have increasingly been used to 
investigate population dynamics on the global scale by means 
of numerical integration of the corresponding system of ODEs 
[iniiiiiiiiE]. Although such investigations captured a rich 
variety of types of behaviour such as fold bifurcations when 
points of equilibrium merge/annihilate [M], or various types of 
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chaotic dynamics |15l I16 |. they provide little analytic insight 
and are limited to small to medium-sized systems. 

In this paper we attempt to investigate the generic prop¬ 
erties of the global dynamics of large complex multi-species 
systems by retaining in our model only the bare essentials - 
nonlinearity and stochasticity. Much in the spirit of May’s 
original approach, the model we propose is simple enough to 
allow for an analytic treatment yet at the same time is rich 
enough to exhibit a non-trivial behaviour. In particular, our 
model captures an instability transition of the May-Wigner 
type, but now on the global scale. It also sheds additional 
light on the nature of this transition by relating it to an ex¬ 
ponential explosion in the number of equilibria. Interestingly, 
despite the nonlinear setting of the problem the random ma¬ 
trix ideas again play a central role in our analysis. 

Similar to the May’s linear model our toy model is likely 
to have rather limited practical significance for quantitative 
description of real ecosystems but it might provide insight into 
the generic qualitative features of such systems and beyond, 
e.g. machine learning [m or financial ecosystems [HIS]. 
The idea of destabilisation by interaction is of relevance far 
beyond the mathematical ecology context as applications of 
systems of many coupled non-linear ODEs are vast (e.g. com¬ 
plex gene regulatory networks [20] , neural networks [211122j , 
or random catalytic reaction networks |23|i. 


Model 

Consider a system of N coupled non-linear autonomous ODEs 
of the form 


-^ =-fiXi + .. ,Xn), i = [2] 

where > 0 and the components /i(x) of the vector field 
f = (/i,..., /jv) are zero mean random functions of the state 
vector X = {xi,... ,xn)- To put this model in the context 
of the discussion above, if Xe is an equilibrium of i.e., if 
+ f(x^ = 0 then, in the immediate neighbourhood of 
Xe system S reduces to May’s model 0 with y = X — Xe 
and Jjk = (ofj/dxk){xe). 

The non-linear system Q may have multiple equilibria 
whose number and locations aepend on the realisation of the 
random field f(x). To visualise the global picture, it is helpful 
to consider first a special case of a gradient-descent flow, char¬ 
acterised by the existence of a potential function V(x) such 
that f = —W. In this case, system ® can be rewritten as 
dx/dt = —\7L, with L{x) — /r|x|^/2 -|-V(x) being the associ¬ 
ated Lyapunov function describing the effective landscape. In 
the domain of L, the state vector x(t) moves in the direction 
of the steepest descent, i.e., perpendicular to the level surfaces 
I/(x) = h towards ever smaller values of h. This provides a 
useful geometric intuition. The term /r|xp/2 represents the 
globally confining parabolic potential, i.e., a deep well on the 
surface of L(x), which does not allow x to escape to infinity. 
At the same time the random potential I^(x) may generate 
many local minima of L{x) (shallow wells) which will play the 
role of attractors for our dynamical system. Moreover, if the 
confining term is strong enough then the full landscape will 
only be a small perturbation of the parabolic well, typically 
with a single stable equilibrium located close to x = 0. In the 
opposite case of relatively weak confining term, the disorder- 
dominated landscape will be characterised by a complicated 
random topology with many points of equilibria, both stable 
and unstable. Note that in physics, complicated energy land¬ 
scapes is a generic feature of glassy systems with intriguingly 


slow long-time relaxation and non-equilibrium dynamics, see 

e.g. [H]. 

The above picture of a gradient-descent flow is however 
only a very special case since the generic systems of ODEs 0 
are not gradient. The latter point can easily be understood in 
the context of model ecosystems. For, by linearising a gradi¬ 
ent flow in a vicinity of any equilibrium, one always obtains a 
symmetric community matrix, whilst the community matrices 
of model ecosystems are in general asymmetric. Note also a 
discussion of an interplay between non-gradient dynamics in 
random environment and glassy behaviour in |25| . 

To allow for a suitable level of generality we therefore sug¬ 
gest to choose the A—dimensional vector field f(x) as a sum 
of ‘gradient’ and non-gradient (‘solenoidal’) contributions: 


/i(x) 


dV{y^) ^ dAjj (x) 

dxi y/N ^ dxj 


i = l,...,A, [3] 


where we require the matrix A(x) to be antisymmetric: Aij = 
—Aji. The meaning of this decomposition is that vector fields 
can be generically divided into a conservative irrotational 
component, sometimes called ‘longitudinal’, whose gradient 
connects the attractors or repellers and a solenoidal curl field, 
also called ‘transversal’. As discussed in, e.g., [55] such a rep¬ 
resentation is closely related to the so-called Hodge decom¬ 
position of differential forms and generalises the well-known 
Helmholtz decomposition of the 3—dimensional vector fields 
into curl-free and divergence-free parts to higher dimensions. 
Correspondingly, we will call V (x) the scalar potential and 
the matrix A(x) the vector potential. The normalising factor 
1 / y/N in front of the sum on the right-hand side in 0 ensures 
that the transversal and longitudinal parts of f(x) are of the 
same order of magnitude for large N. 

Finally, to make the model as simple as possible and 
amenable to a rigorous and detailed mathematical analysis we 
choose the scalar potential V (x) and the components Aij (x), 
i < j, of the vector potential to be statistically independent, 
zero mean Gaussian random fields, with smooth realisations 
and the additional assumptions of homogeneity (translational 
invariance) and isotropy reflected in the covariance structure: 

(y(x)l/(y)) - i;^rv(|x-y|^) ; [4] 

(x)-AriTTT, (y)) — d r^l^lx y| ^{^SinSjrn • [ 5 ] 

Here the angular brackets (...) stand for the ensemble average 
over all realisations of V (x) and A(x), and Sin is the Kronecker 
delta: Sin = 1 if f = n and zero otherwise. 

For simplicity, we also assume that the functions rv(r) 
and Fa)?") do not depend on N. This implies [27] 

POO 

ra{r) = / exp(—sr) 7 ,j(s)ds, a ^ A,V, 

Jo 

where the ‘radial spectral’ densities 7cr(s) > 0 have finite to¬ 
tal mass: 7 o-(s)ds < oo. We normalize these densities by 

requiring that r"(0) = s^ 7 cr(s)ds = 1. The ratio 

T = u^/(u^ + a^), 0 < r < 1, 

is a dimensionless measure of the relative strengths of the lon¬ 
gitudinal and transversal components of f(x): if r = 0 then 
f(x) is divergence free and if r = 1 it is curl free. 


Results 

Determining and classifying all points of equilibria of a dy¬ 
namical system with many degrees of freedom is a well-known 
formidable analytical and computational problem. In this pa¬ 
per we shall focus our investigation on the simplest, yet in¬ 
formative characteristic of system 0 by counting its total 
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number of equilibria, that is the total number Ntot of solu¬ 
tions of the simultaneous equations 

- -I- .. .,xn) = 0, i = I,... ,N. [6] 

Certainly, finding Aftot is a good starting point of any phase 
portrait analysis. 

Had we restricted ourselves to the gradient-descent flows, 
Aftot would simply count the number of stationary points 
(minima, maxima, or saddle-points) on the surface of the Lya¬ 
punov function L(x). The problem of counting and classifying 
stationary points of high-dimensional random energy land¬ 
scapes of various types attracted considerable interest in re¬ 
cent years In particular, works |28l I29| 

study such energy landscapes generated by a potential equiv¬ 
alent to the above Lyapunov function. One of the main con¬ 
clusions of that study is that for N large the topology of the 
Lyapunov function changes drastically with decrease of the 
strength of the confining term relative to that of the interac¬ 
tion term in L(x). The change manifests itself in the emer¬ 
gence of multitude of equilibria, exponential in number. Such 
a transition is intimately connected to the spin-glass like re¬ 
structuring of the Boltzmann-Gibbs measure induced by the 
Lyapunov function when the latter is treated as an effective 
energy landscape. 

We shall prove below that for N large the general au¬ 
tonomous system S exhibits a similar drastic change in 
the total number of equilibria when the control parameter 

m = —where a = 2 \/+ a'^, 
aVN 

drops below the threshold value rric = 1. As in the case of 
gradient systems, the proof involves the Kac-Rice formula as 
a starting point. However, performing the subsequent steps 
requires quite different mathematical techniques due to the 
asymmetry of the Jacobian matrix for non-gradient systems. 

The Kac-Rice formula, see e.g. [34], counts solutions of si¬ 
multaneous algebraic equations. Under our assumptions (ho¬ 
mogeneity, isotropy and Gaussianity of V and A ), this formula 
yields the ensemble average of Ntot in terms of that of the 
modulus of the spectral determinant of the Jacobian matrix 
Jij = dfijdxj (see Materials and Methods): 

(Aftot) = ^ (|det(-^<5ij-h Jij)!), [7] 

A* 

thus bringing the original non-linear problem into the realms 
of the random matrix theory. 

The probability (ensemble) distribution of the matrix J 
can easily be determined in closed form. Indeed, the matrix 
entries of J are zero mean Gaussian variables and their co- 
variance structure, at spatial point x, can be obtained from 
0-0 by differentiation: 

{Jijtlnm) — Oi [(1 -f (x jn^irn J- JijJmn)] , 

where ejv = (1 — t)/N. Thus, to leading order in the limit 
N oo, 

Jij = a{Xij + y/rSiji), [ 8 ] 

where Xij^ i,j = 1,..., N are zero mean Gaussians with 

i^Xij Xnrn) — ^in^jm + T5jnSim, [ 9 ] 

and ^ is a standard Gaussian, ^ ~ N{G, 1), which is statisti¬ 
cally independent of X = {Xtj). Note that for the divergence 
free fields f(x) (i.e., if r = 0) the entries of J are statisti¬ 
cally independent in the limit X —>■ oo, exactly as in May’s 


model. On the other side, if f(x) has a longitudinal compo¬ 
nent (r > 0) then this implies positive correlation between 
the pairs of matrix entries of J symmetric about the main 
diagonal: (XijXji) = t ii i j. Such distributions of the 
community matrix has also been used in the neighbourhood 
stability analysis of model ecosystems [S]. Finally, in the lim¬ 
iting case of curl free fields (r = 1), the matrix J is real 
symmetric. 

The representation ^ comes in handy as it allows one to 
express 0 as a random matrix integral: 


{Not) — 


N- 


f 


_ -N't , 

e ^ dtt 


(101 


I det (x 5ij Xij ) I) ■^ 2-kIN ' 

where x = \fN{m J- t\/r) and the angle brackets (.. .)xjv 
stand for averaging over the real elliptic ense mbl e of random 
N X N matrices X defined in §, see also ( |20[ ). This one- 
parameter family of random matrices interpolates between 
the Gaussian Orthogonal Ensemble of real symmetric matrices 
(GOE, r = 1) and real Ginibre ensemble of fully asymmetric 
matrices (rGinE, r = 0), see [35] for discussions. Both rGinE 
and its one-parameter extension 0 have enjoyed considerable 
interest in the literature in recent years [SllSTiiiiiiiiiinj. 

The matrix X is asymmetric (unless r = 1) and can 
have real as well as complex eigenvalues. The latter come in 
complex-conjugate pairs. Their density, in the limit N —>■ oo, 
is constant inside the ellipse with the main half-axis NN (1±t) 
and vanishes sharply outside [4111391135] . The corresponding 
theorem is known as the Elliptic Law and its validity extends 
beyond the Gaussian matrix distributions [421143] . However, 
in the context of our investigation it is the density of real 
eigenvalues of X that appears to be most relevant. 

Denote by (x) the density of real eigenvalues ot N x N 
matrices X ® averaged over all realisations of X. It is con¬ 
venient to normalize ff^\x) in such a way that Pn{x) dx 
gives the average number of real eigenvalues of X in the in¬ 
terval [q:,/3]. a crucial observation is that p^^(x) is directly 
related to the averag ed v alue of the modulus of the determi¬ 
nant that appears in ( jlOj ). Namely, 

(|det(xJij-Xij)|)_y^ =Civ(r)e2(f+-) p%l^{x), [11] 

where Cn{t) = 2^/l + t {N — 1)!/(X — 2)!! and p^^]^(x) is 
the average density of real eigenvalues of matrices X of size 
{N J- 1) X (N + 1). For the limiting case r = 0 this rela¬ 
tion appeared originally in )44| . and it can be extended to any 
r £ [0,1) without much difficulty (see SI for a derivation of 
( jllj ) following the approach of m )• In the limiting case of 
real symmetr ic m atrices r = 1, all eigenvalues of X are real 
and relation dll]) is also vali d ]28 1. 

Combining |lo| and (111 and changing the variable of in¬ 
tegration from t to X = m + t^/r, one can express {Ntot) for 
system 0 with N degrees of freedom in terms of the density 
of real eigenvalues in the elliptic ensemble of random matrices 
0 of size (X J- 1) X {N J- 1): 


{Not) — 


1Cn{t) 


f 


e--®(^)piv+i(AVx)^, [12] 

o V 


where S(A) = 2 ^^ ~ 2 {f+T) andXjv(r) = N 2 ^ CN{T)/y/T. 

The importance of this relation is due to the fact that p^^ (x) 
is known in closed form in terms of Hermite polynomials [ 39] . 
This allow s us to carry out an asymptotic evaluation of the 
integral in (121 and calculate {Not) in the limit N —>■ 00 . The 
key finding that emerges from this calculation is that {Not) 
changes drastically around m = 1. If m > 1 then 


lim {Not) = 1. 

JV^-oo 


[13] 
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On the other hand, if 0 < m < 1 then, to leading order in the 
limit —>■ 00 , 


{Ntot) = 


fl4l 


0 < u < oo, one obtains 7 r = f^- 


where Etot(m) = |(jn^ — 1) — In m > 0 for all 0 < m < 1. 
Therefore, if m < 1 then {Ntot) grows exponentially with 
N. The factor in front of the exponential in (141 is given by 
7 t = V2(l + r)/(l-r) as long as r < 1. The gradient limit 
T = 1 can be approached by scaling r with N . Setting t = 1 — 

e-'^^p^dp. 

This regime describes a weakly non-gradient flow. The cor¬ 
responding regime for ensembles of asymmetric matrices was 
discovered long ago uni En¬ 
close to the phase transition point m = 1 the complexity 
exponent vanishes quadratically, Etot = (1 — m)^ as m —>■ 1, 
implying that the width of the transition region around m — 1 
is l/\/iV. According to the general lore of phase transi¬ 
tions, for large but finite N there exists a ‘critical regime’ 
m = 1 + where the num ber of e quili bria changes 

smoothly b etwe en the two ‘phases’ (131 and (141. A quick in¬ 
spection of (121 shows that the corresponding crossover profile 
is determined by the profile of in the vicinity of the 


‘spectral edge’ a; = (1 -|- t)-\/N, see Materials and Methods. 
After rescaling A, A = 1 -|-t-|- , the density p^^(A\/iV) 


converges to , ^ pidleCC) in the limit N ^ oo, where [39]: 

\/ 1 — T^ ^ 


pid 9 e(C) = |erfc(V2C) + ‘i" [1 -f erf (C)]| [15] 

with erf (x) = 1—erfc(a;) = fo In terms of 

the critical crossover profile is given by 

(Mot) = 7r e"T pM ^ + k dt [16] 

where Ct = \/t/(1 — r). The right-hand side of (|l 6|) i nterpo- 
lates smoothly between the two regimes ( |13| ) and ( |i4| |, when 
parameter k runs from k = —oo to k = -|-oo. 

Although our investigation is concerned with the ensem¬ 
ble average of the number of equilibria Mtot, we expect that 
in the limit N ^ oo the deviations of Mtot from its average 
(Mtot) are relatively small. This is certainly the case above 
the critical threshold, for m > 1. For, under some additional 
technical assumptions on the decay of correlations for f(x), 
the system © will almost certainly have at least one station¬ 
ary point, see ES] for the relevant results about the maxima 
of homogeneous Gaussian fields. Therefore, Mtot > 1 and the 
established convergence of (Mtot) to 1 in the limit N ^ oo 
actually implies that the probability for Mtot to take other 
values than one is close to zero for large N. The problem of 
estimating the deviation of Mtot from its average value in the 
opposite regime 0 < m < 1 is much harder and is an open 
and interesting questiorj^ 


Discussion 

In this paper we introduced a model describing generic large 
complex systems and examined the dependence of the total 
number of equilibria in such systems on the system complex¬ 
ity as measured by the number of degrees of freedom and the 
interaction strength. The inspiration for our work came from 
May’s pioneering study [T] of the neighbourhood stability of 
large model ecosystems. Our outlook is complementary to 
that of May’s in that it adopts a global point of view which 


is not limited to the neighbourhood of the presumed equilib- 


In the context of model ecosystems our analysis is applica¬ 
ble to complex multi-species communities in which each kind 
of species on its own becomes extinct and thus interaction be¬ 
tween species is key to persistence of the community. The key 
feature of our analysis is that in the presence of interactions, 
as the complexity increases, there is an abrupt change from 
a simple set of equilibria (typically, a single equilibrium for 
large number of species ^ 1) to a complex set of equilib¬ 
ria, with their total number growing exponentially with N. In 
the latter regime, we expect the stable equilibria to be only a 
tiny proportion of all the multitude of equilibria, see discus¬ 
sion below, which is indicative of long relaxation times and 
transient non-equilibrium behaviour. 

We expect this sharp transition in the phase portrait to 
be shared by other systems of randomly coupled autonomous 
ode’s with large number of degrees of freedom. To that end, 
it is appropriate to mention that very recently a similar ‘explo¬ 
sion in complexity’ was reported in a model of neural network 
consisting of randomly interconnected neural units ^22] ■ The 
model considered in [22] is essentially of the form (|^ but with 
the particular choice of ft = JijS{xj) where o is an odd 
sigmoid function representing the synaptic nonlinearity and 
Jij are independent centred Gaussian variables representing 
the synaptic connectivity between neuron i and j. Although 
being Gaussian, the corresponding (non-gradient) vector field 
is not homogeneous and thus seems rather different from our 
choice and not easily amenable to a rigorous analysis. Never¬ 
theless, a shrewd semi-heuristic analysis of [2^ revealed that 
close to the critical coupling threshold the two models actu¬ 
ally display very similar behaviour, described essentially by 
the same exponential growth in the total number of equilibria 
with rate Etot(m). This fact po ints tow ard s considerable uni¬ 
versality of the transiti on f rom (131 to (141 and suggests that 
the crossover function (161 may be universal as well. 

Our model captures an abrupt change in the dynamics 
of large complex systems on the macroscopic scale. At the 
same time zooming in to classify each and every equilibrium 
point into locally stable or unstable seems a hard task. For, 
although linearising the field f(x) around a given equilibrium 
is fairly straightforward, with the outcome being the Jacobian 
matrix ([^, conditioning by the positions of equilibria and tak¬ 
ing into account all eventualities is a highly non-trivial task. 
Given the stochastic setup of our model the question about 
stability of individual equilibria may be even the wrong ques¬ 
tion to ask, whereas addressing the statistics of the number of 
stable equilibria seems very appropriate. 

Arguments similar to those in the previous section yield 
the ensemble average of the total number of stable equilibria, 
(Mst), over all realisations o f th e vector field f(x) in terms the 
random matrix integral (cf.lUM): 


(Mst) = 


N- 


OO 

J {det(a: 


Sij Xij)xxiXij)) 




_ Nt^ , 

e 2 dt 

yWA ’ 


[17] 


where Xx(X) = 1, if all A eigenvalues of matrix X have real 
parts less than the spectral parameter x = MN{m + t^/T)j 
XxiX) = 0 otherwise. In the limiting case of a purely gradi¬ 
ent dynamics t = 1, the rescaled Jacobian matrix X is real 


-^In this context we would like to mention the recent work of Subag 1331 who proved that the 
deviations of A/tot from (A/tot) in the spherical p-spin model are negligible in the limit of large 
system size. Though that model is different from ours, it is not dissimilar to the gradient limit 
of T = 1 of our model [3^, for instance, the average number of equilibria grows exponentially 
with N 1301. T hus one might hope to adopt the technique of 1331 to our model. Another relevant 
reference is 1491 . 
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symmetric with all N eigenvalues real. In this case the above 
integral can be related to the probability density of the max¬ 
imal eigenvalue of the GOE matrix |29l [30], with the latter 
being a well-studied object in the random matrix theory, see 
e.g. and references therein. This observation can then 

be used to evaluate (Afst) for N ^ 1. One finds [55] that 
(Afst) —>■ 1 if m > 1, whilst if 0 < m < 1 then, to leading 
order in N, (Afst) oc with 0 < Tist < Thus, in the 

case of purely gradient dynamics, as the complexity increases, 
large nonlinear autonomous systems assembled at random un¬ 
dergo an abrupt change from a typical phase portrait with a 
single stable equilibrium to a phase portrait dominated by 
an exponential number of unstable equilibria with an admix¬ 
ture of a smaller, but still exponential in N, number of stable 
equilibria. 

It was suggested to us by J.-P. Bouchaud that in the gen¬ 
eral case of non-gradient dynamics 0 < t < 1, it would be nat¬ 
ural to expect a further phase transition in the plane (m, r) 
such that below a certain number rc{m) stable equilibria are 
no longer exponentially abundant in the limit N ^ oo (i.e. 
'Sst(m,T) —>■ 0), with further implications for the global dy¬ 
namics. Unfortunately, for a fixed 0 < r < 1 only vanishing 


fraction of order of of eigenvalue s of X remain real, 

and the relation of the integral in Eq. (171 to statistics of 
the largest real eigenvalue in the elliptic ensemble seems to 
be lost. This fact prevented us so far from reliable counting 
of stable equilibria for the general case of non-gradient flows. 
In principle, for given values of parameters N, r, m one may 
attempt to evaluate the ensemble average in the integral in 
Eq. ( |17[ ) numerically, and then to evaluate numerically the 
integral itself. Although such a procedure seems tractable, its 
actual implementation with sufficient precision is not straight¬ 
forward, especially in the limit N ^ oo due to the exponen¬ 
tially large values involved. Clarification of the status of the 
picture suggested by J.-P. Bouchaud and related studies re¬ 
main an important outstanding issue and is left for a future 
work. 


Materials and Methods 

Kac-Rice Formula The expected number {Aftot) of simultaneous solutions to the 
system of equations in is given by the formula, see e.g. 1341 . 

{Mtot)= f ((5(-/rx+f(x)) |det (-p<5ij-|-Jij(x))|) , [18] 

Jr" 

where 5(x) is the multivariate Dirac (5-function, dx^ is the volume element in 
and Jij[x) = dfijdxj are matrix elements of the Jacobian matrix J = (t/ij). 
By our assumptions - ||^ the random field f(x) is homogeneous and isotropic. 
For such fields samples of I and J taken in one and the same spatial point X are 
uncorrelated, {fi • dfijdXj) = 0 for all This is well known and can be 

checked by straightforward differentiation. In addition, the field f is Gaussian, hence 
the f(x) and J(x) are actually statistically independent. This simplifies the eval¬ 
uation of the integral in | |18| considerably. Indeed, the statistical average in jl8| 
factorizes and, and since (| det {—fiSij + Jij{x)) |) does not vary with X, the 
integrand effectively reduces to 

__^-^(k.x(^,k.f(x))^ [ 19 ] 

Furthermore, at every spatial point X the vector f(x) is Gaussian with uncorrelated 
and identically distributed components, 

(A (x)/, (x)> = , (T^ = 2v^ irV (0) I -F 2a2 (0) I . 

Therefore , and evaluating the integral on the right-hand 

side in \19) one arrives at Q. 


random matrices X of size X is given by 

1 


Pjv(A) = .2-1 exp [- 


2(1-r2) 


Tr {XX^ - tX‘ 


[ 20 ] 


where 2^ is the normalisation constant and r £ [0, 1). It is straightforward to 
verify that the covariance of matrix entries is given by the expression specified 

i n ^ 9). The mean density of real eigenvalues of in the elliptic ensemble 

polls known in closed form in terms of Hermite polynomials, see 1391 . Assuming 
for simplicity that A^ + 1 is even, one has 
where 

N-1 I , (t) / u 2 

_ _L V [21] 




k=0 


k\ 


and 


\/2^’(l-|-r)(Af-1)! 


f [ 22 ] 

Jo 


Here (x) = e {x) and h^jj\x) are rescaled Hermite polyno¬ 
mials, h^\x) = -^= dt. This, together with the 

integral | |l2| allow one to evaluate {J\f)tot in the limit N —>■ CXD. We shall sketch 
the corresponding evaluation below. 

The asymptotics of p^^ (a;) in the bulk and at the edge of the support of the 
distribution of real eigenvalues in the real elliptic ensemble were found in 1391 , and 


outside the support it can also be readily extracted using HD 
in the bulk, i.e., for |a:| < (1 + the contribution of 

dominant and, to leading order in N, 

1 


P^ii(AVA)| 


22^ . In particular, 
^ to p%\x) is 

[23] 


y'27r(l - t2) ■ 

At the same time for |a;| > (1 + T)\/iV both Ih) and [ 22 ] yield exponentially 
small contributions to with |22| being dominant. Our evaluation yields 

^ = Q(A)exp-Af'J(A), where 


Q(A) = 


27r(l + t) VA2 - 4r)(A - y/X^ - 4r) 


1/2 


[24] 


A^ 


[25] 


= 2{rTT ^^ 2 ^ "" 

The form of HD suggest the application of the Laplace method. One easily finds that 
5'{ A) has a minimum at A* = m(lH-T) which belongs to the domain | A| < 1 H-t 
as long as 0 < 771 < l.Thus, applying the Laplace method and taking into account 
the asymptotic formula 

/Civ(r) + [26] 

one arrives at the asymptotic expression jl4| for {J^)tot in ths parameter range 
0 < m < 1. For m > 1 the saddle-point occurs in the domain A > 1 + r so 
that the analysis requires search for the minimum of S'(A) + 'I^{A). After straigh- 
forward algebra we find ^ (^{A) + ^'(A)) = ~4 t) — ^ which 

is equal to zero atA = A* = m ^ One also verifies that this is 

a point of minimum for ^(A) + ^(A) and a further simple calculation then yields 
^{A*) + ^(A*) = — ln(m/\/r). Calculating the saddle-point contribution 
then yields jl3[ |. 

The above asymptotic analysis assumes that 0 < t < 1. Let us now discuss 
the modifications required to study the scaling regime of weakly non-gradient flow 
T —¥ 1 for 0 < TTL < 1. We only need to evaluate the leading contribution to 
p^_l_^ given by |21|. By making use of the above integral representation for the 

scaled Hermite polynomials (x) and applying the scaling t = 1 — ^, we can 
write 

with 7jv(A) given by 


Real elliptic matrices and asymptotics of {J\f)tot The joint probability den¬ 
sity function of the matrix entries in the elliptic ensemble of real Gaussian 






^N-2 
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where $ 7 ^( 0 .) = e Recalling that the limit of $;v(a) as 

A^^ooislifO<a<l and 0 if a > 1, one obtains 


P^Vi(aV^) 



dp. 


A = 1 + r 


Cv^ 


+ T H-—’ where we have, to the leading order in N, 


„-iV(S(A)+i) 


= exp 


1 — T 
2 r 


(^" + C') + 


VT^- 


-kQ 


This together with j26^ and (iD converts IH) to 


Substituting this expression into the integrand of IH) and evaluating the integral in 
the limit N ^ 00 (hence, r ^ 1) by the Laplace method then yields {J\f)tot 
the weakly non-gradient regime. 

Finally, our calculation of the profile of {J\f)tot in the transitional region 
m = 1 + uses the fact that in such a regime the main contri¬ 

bution to the integral (HI comes from the neighbourhood of the spectral edge, 


ACKNOWLEDGMENTS. The research of the first author was supported by EPSRC 
grant EP/J002763/1 “Insights into Disordered Landscapes via Random Matrix Theory 
and Statistical Mechanics". The second author would like to thank the Isaac Newton 
Institute for Mathematical Sciences for its hospitality during the programme Periodic 
and Ergodic Spectral Problems supported by EPSRC Grant Number EP/K032208/1. 
Both authors thank David Arrowsmith, Matthew Evans and Vito Latora for com¬ 
ments on various versions of this paper and to J.-P. Bouchaud for his comments on 
the number of stable equilibria in non-gradient systems. 


1. May RM (1972) Will a Large Complex System be Stable? Nature 238: 413-414. 

2. Gardner MR, Ashby WR (1970) Connectance of Large Dynamic (Cybernetic) Systems: 
Critical Values for Stability Nature 228: 784. 

3. Porter CE (1965) Statistical Theories of Spectra: Fluctuations. A Collection of 
Reprints and Original Papers (Academic Press, New York). 

4. Allesina S, Tang Si (2015) The stability complexity relationship at age 40: a random 
matrix perspective. Popul Ecol 57: 63-75. 

5. James A, Plank MJ, Rossberg AG, Beecham J, Emmerson M, and Pitchford JW 
(2015) Constructing Random Matrices to Represent Real Ecosystems. Am Nat 185: 
680-692. 

6. McCann KS (2000) The diversity-stability debate. Nature 405: 228-233. 

7. Boyd (2012) The art of ecological modeling. Science 337: 306-307. 

8. Allesina S, Tang Si (2012) Stability criteria for complex ecosystems. Nature 483: 205 
-208. 

9. Johnson S, Domnguez-Garca V, Donetti L, Muoz MA (2014) Trophic coherence de¬ 
termines food-web stability. Proc Natl Acad Sci USA 111: 17923-17928. 

10. Yodzis P, Innes S (1992) Body-size and consumer-resource dynamics. Am Nat 139: 
1151-1173. 

11. Williams RJ, Martinez ND (2004) Stabilization of chaotic and non-permanent food- 
web dynamics. European Physical Journal B 38: 297-303. 

12. Dunne JA, Brose U, Williams RJ. Martinez ND (2005) Modeling food-web dynam¬ 
ics: complexity-stability implications. In Aquatic Food Webs: An ecosystem approach, 
OUP 2005. 

13. Rohr RP, Saavedra S, Bascompte J (2014) On the structural stability of mutualistic 
systems. Science 345: 1253497. 

14. Dai L, Vorselen D. Korolev KS. Gore J (2012) Generic indicators for loss of resilience 
before a tipping point leading to population collapse. Science 336: 1175. 

15. Hastings A, Powell T (1991) Chaos in a three-species food chain. Ecology 72: 896- 
903. 

16. McCann K, Yodzis P (1994) Non-linear dynamics and population disappearances . Am 
Nat 144: 873-879. 

17. Galla T, Farmer JD (2013) Complex dynamics in learning complicated games. Proc 
Natl Acad Sci USA 110: 1232-1236. 

18. Haldane AG. May RM (2011) Systemic risk in banking ecosystems. Nature 469: 351- 
355. 

19. Farmer JD. Skouras S (2013) An ecological perspective on the future of computer 
trading. Quantitative Finance 13: 325-346. 

20. de Jong H (2002) Modeling and Simulation of Genetic Regulatory Systems: A Liter¬ 
ature Review. J Comput Biol 9: 67-107. 

21. Sompolinsky H, Crisanti A, Sommers H-J (1988) Chaos in Random Neural Networks. 
Phys Rev Lett 61: 259 - 262. 

22. Wainrib G, Touboul J (2013) Topological and Dynamical Complexity of Random Neu¬ 
ral Networks. Phys Rev Lett 110: 118101. 

23. Stadler PF, Fontana W. Miller JH (1993) Random Catalytic Reaction Networks . 
Physica D 63: 378 - 392. 

24. Angelani L., Paris! G., Ruocco G., Viliani G. (2000) Potential Energy Landscape and 
long-time dynamics in a simple model glass. Phys Rev E 61: 1681 - 1691. 

25. Cugliandolo LF , Kurchan J, Le Doussal P, and Peliti L (1997) Glassy behaviour in 
disordered systems with nonrelaxational dynamics Phys Rev Lett 78: 350-353. 

26. Zhou J Xu, Allyu MDS, Aurell E, Huang Sui (2012) Quasi-potential Landscape in 
complex multi-stable systems J R Soc Interface 9: 3539 - 3553. 

27. Schoenberg IG (1938) Metric Spaces and Completely Monotone Functions. Ann Math 
39: 811-841. 

28. Fyodorov YV (2004) Complexity of Random Energy Landscapes, Glass Transition, and 
Absolute Value of the Spectral Determinant of Random Matrices. Phys Rev Lett 92: 
240601. 

29. Fyodorov YV, Nadal C (2012) Critical Behavior of the Number of Minima of a Random 
Landscape at the Glass Transition Point and the Tracy-Widom Distribution. Phys Rev 
Lett 109: 167203. 

30. Aufhnger A, Ben Arous G, Cerny J (2013) Random Matrices and Complexity of Spin 
Glasses. Common Pure AppI Math 66(2): 165 - 201; 


31. Nicolaescu L (2014) Complexity of random smooth functions on compact manifolds. 
Indiana University Mathematics Journal 63(4): 1037-1065 

32. Fyodorov YV (2015) High-Dimensional Random Fields and Random Matrix Theory. 
Markov Processes and Related Fields 21(3): 483-518. 

33. Subag E (2015) The complexity of spherical p-spin models - a second moment ap¬ 
proach. arXiv:1504.0225 

34. Azais J-M, Wschebor M (2009) Level Sets and Extrema of Random Processes and 
Fields (John Wiley &£. Sons). 

35. Khoruzhenko BA and Sommers H-J (2011) Non-Hermitian Ensembles. Chapter 18 of 
The Oxford Handbook of Random Matrix Theory, eds G. Akemann, J. Baik and P. Di 
Francesco (Oxford University Press), pp 376-397. 

36. Akemann G, Kanzieper E (2007) Integrable structure of Ginibre’s ensemble of random 
real matrices. J Stat Phys 129: 1159-1231. 

37. Sommers HJ, Wieczorek W. (2008) General eigenvalue corrrelations for the real Gini- 
bre ensemble. J Phys A Math Theor 41: 405003. 

38. Borodin A, Sinclair CD (2009) The Ginibre Ensemble of Real Random Matrices and 
its Scaling Limits. Communications in Mathematical Physics 291: 177-124. 

39. Forrester P J, Nagao T (2008) Skew orthogonal polynomials and the partly symmetric 
real Ginibre ensemble. J Phys A Math Theor 41: 375003. 

40. Akemann G, Phillips MJ (2014) The Interpolating Airy Kernels for the = l and 
0 = 4 Elliptic Ginibre Ensembles. J Stat Phys 155: 421 - 465. 

41. Sommers H-J, Crisanty A, Sompolinsky H, Stein Y (1998). Spectrum of large random 
asymmetric matrices Phys Rev Lett 60: 1895-1898. 

42. Girko VL (1995) The Elliptic Law: Ten years later Random Operators and Stochastic 
Equations 3(3): 257-302. 

43. Nguyen H and O’Rourke S (2015) The elliptic law. International Mathematics Re¬ 
search Notices 2015(17): 7620-7689. 

44. Edelman A, Kostlan E, Shub M (1994) How many eigenvalues of a real matrix are 
real? Journal of the American Mathematical Society 7: 247 - 267. 

45. Fyodorov YV and Khoruzhenko BA (2007) On absolute moments of characteristic 
polynomials of a certain class of complex random matrices. Communications in Math¬ 
ematical Physics 273: 561 - 599. 

46. Fyodorov YV, Khoruzhenko BA, Sommers H-J (1997) Almost-Hermitian random ma¬ 
trices: Eigenvalue density in the complex plane. Phys Lett A 226: 46 - 52. 

47. Efetov KB (1997) Directed quantum chaos. Phys Rev Lett 79: 491 - 494. 

48. Piterbarg, VI (1996) Asymptotic Methods in the Theory of Gaussian Processes and 
Fields (American Mathematical Society, Providence). 

49. Nicolaescu L (2015) A CLT concerning critical points of random functions on a Eu¬ 
clidean space. arXiv:1509.06200 

50. Majumdar SN, Schehr G (2014) Top eigenvalue of a random matrix: large deviations 
and third order phase transition. J Stat Mech P01012. 


Appendix: Supporting Information 

In this appendix we express the mean density o f rea l eigen¬ 
values in the real elliptic ensemble, see Equation (351 below, 
of X matrices Xn in terms of the modulus of the spectral 
determinant of real elliptic matrices of size {N — 1) x {N — 1). 
Our derivation is based on the method suggested in [44], how¬ 
ever our Jacobian computation differs from the one given in 
[44] and may be of interest on its own. For a similar calcula¬ 
tion in the context of complex matrices see [45]. 
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Housholder reflections and partial triangularization of real 
matrices The key idea of [44] is based on employing House¬ 
holder reflections described by matrices 


Pv = /jv — 2v ( 


T 1 

V V = 1. 


[271 


where In is N x N identity matrix, v is column vector in 
of unit length. Its transpose, v^, is row vector, so that 
the Kronecker product v (g) is a matrix. Pv describes the 
reflection about the hyperplane with normal v and passing 
through the origin. Obviously, Pv is is symmetric and orthog¬ 
onal: Pj = Pv and Pv = In- 

Let ei be the first vector of the standard basis in R^, i.e., 
ef = (1,0,..., 0). For any unit vector x, jxj = 1 define 


X -f ei 

V'2(1+2;i)' 


[281 


and consider the Hausholder reflection Pv built from the above 


V according to ( 271 . Then it is easy to check that PvX = —ei. 


This implies that for any nonzero vector x there exists a 
Housholder reflection such that PvX = fcei, with k = —jxj. 

Let A be a real eigenvalue of the N x N matrix with 
real entries A[^\ i.e. = Ax for some column A'’—vector 

X of unit length. Our goal is to demonstrate that it is always 
possible to represent that matrix as 


4(Jv)_pA 


[291 


for some real {N — 1)- component vector w and a real matrix 
X (jY —1). To verify this take the Housh¬ 
older reflection acting on the eigenvector x as PvX = fcei. Ob¬ 
viously PvA’-^^x = APvX = Afeei. In view of P^ = In we can 
rewrite the left-hand side as P^A^^^PvPvX = k PvA^’^^P^ei 
which after denoting A = PvA^^^Pv imples Aei = Aei. Us¬ 
ing the definition of ei we see that An = A and Aij = 0 
for all j = 2,..., N. Therefore A can be written as ^4 = 

and as A^^'i = PvAPv the relation (291 fol- 


0 

lows. 


Considering now the volume element dA^^^ = dA^^'' 

our next goal to write it down in terms of in dependent vari¬ 


ables parametrizing the right-hand side of (291, that is (A'^—1)^ 
variables parametrizing A'^ — 1 components of w, 1 

parameter for A, and the remaining A'^ — 1 parameters for rep¬ 
resenting the matrix P . A convenient parametrization for P 
comes from employing ( |28| l for the vector v, which shows that 
the last A^ — 1 components of that vector {v 2 , ■ ■ ■ ,vnY' = 9 
can be used as independent variables, whereas normalization 

fixes the first component. Writing = ^-y/l — q^q, with 
jqj < 1 and employing (271 yields an explicit parametrization 


P = 


-2q^V^l - q^q 
-2q\/r^-q^ liv-i - 2q(g)q^ 


2q^q - 1 


[30] 


The problem therefore amounts to calculating the Jacobian of 
the transformation » (A , w, q, A^^ ^)). To that end we 

start with differentiating (|29|l which gives 


dA^’^^=P 


(PdP), 


0 


W 


dX 

0 




where we employed that dP P = —P dP and used the nota¬ 
tion [A, E\ = AB — B A for the matrix commutator. A direct 
calculation using (|30|l shows that the matrix (PdP) can be 


symbolically written as 


(PdP) = 


0 

db 


-db^ 

dP 


[31] 


where the expression for dP is immaterial for our goals, and 

=q(Xlq‘ I rfq [32] 


db = ( -2^/1-q^q+ ^ 2q^O„o„T 


\/i - q^q 


-I- f —4\/l — q^q -I- 


1 - 2qS 




q^ q 


(dq^q) q. 


Substituting (311 into the expression for dA^^'i we find that 

,,{N)_-p( dA — w^db dw^-|-db^(A — dP 

'v(A-A(^-i))db dA(^-i)+db(g)w^ + [dP,A(^-i)] 

From this expression we easily read off the required Jacobian 
to be given by 


Jacobian = [ det (A7jv-i — A^^ '^^)[ 


9q 


[33] 


where the last factor symbolically denotes the part of the Ja¬ 
cobian coming from the transformation db —>■ dq described in 
(321. A straightforward calculation shows that 


9b 

9q 


= 2^(l-q^q)^ 


so that finally we arrive at the change-of-measure formula 
dA^"^) = [34] 

2^(1—q^q)T“'^ I det(A7Ar-i—A*-^~^^)[ d\dw^~^dq^~^dA^^~^'^ . 

Elliptic Ensemble of Gaussian Random Matrices. The Joint 
Probability Density (JPD) of the Elliptic Ensemble of Gaus¬ 
sian random matrices Xn of size N x N whose entries have 
covariance (XijXnm) = Si„Sjm + TSjnSim is given by 

T{Xn) = Zn~^ (Tr(AivX^) - r Tr(Alr)) , 

[35] 


'2(l-r2) 


where Zn is the corresponding normalization factor 


Our goal is to find the JPD of variables A, w, q, X jy-i used 
to perform the partial triangulation of Ajv via ( |29| ), with the 
role of A^^^ played now by Xn- We have 

Tr = A'^+w^w-fTr Aiv-iA^_i, Tr A|, = A^'-fTr A^_i . 


iV(iV-l)/4 


Taking into account the change-of-measure formula (341 we 
see that the corresponding JPD can be written as 


P(A,w,q, Aiv-i) = 2^(1 - q'^q) 2 


_1 -T—W W 

e X 


—[det(A7iv-i - Aiv-i)[P(Ajv-i). 

By definition, the density of real eigenvalues p^^(A) is ob¬ 
tained by integrating the above JPD over variables [q[ < 
1,-00 < w < 00 and finally over Av-i. After performing 
the integrals over q and w we immediately arrive at the rela¬ 
tion 

= Cn\{t) ~ . 

where Civ-i(T) is a certain constant which can be found ex¬ 
plicitly. This is precisely equivalent to Equation (13) in the 
main text with the obvious change of notations: A —>■ a: and 
A A-h 1. 
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